
## R code

library(foreign)
library(arm)
library(texreg)
library(ebal)
library(xtable)
library(rjags)
library(ggplot2)
library(mgcv)
library(cem)
set.seed(pi/2)

# detach these packages for now-
# this would cause problems for texreg
detach("package:cem")
detach("package:mgcv")
detach("package:nlme")

###################################
#### Create the final dataset  ###
##################################

#setwd("[[[[REPLICATION ARCHIVE LOCATION HERE]]]")
load(file="cities_140430.RData")

###################################
### Basic Codebook				###
###################################

#######
### Policy variables of interest 
#######

## policy = scaled measure of policy outcomes based on ICMA sustainability dataset

## The following variables are all based on Census of Government 
## salestax_share = Share of taxes from sales tax from Census of Government
## expenditures_capita = Total expenditures per capita from Census of Government
## taxes_capita = Total taxes per capita from Census of Government

#######
### Institutional variables of interest 
#######
## initiative2 -- Factor for whether city has direct democracy (Primarily based on ICMA FoG surveys; cities that were missing in the ICMA FoG files were filled in based on data from the Initiative & Referendum Institute)
## partisan_elections2 -- Factor for whether city has partisan elections (Primarily based on ICMA FoG surveys; cities that were mmissing in the ICMA FoG files were filled in based on CW Internet research)
## offcycle_elections -- Factor for whether city has off cycle elections (based on CW Internet research; by early May there should be more data for this variable)
## council_salary -- City council salaries (from ICMA surveys)
## fog2 -- Factor for whether a city has a mayor-council system or council-manager system (1=mayor)
## perc_atlarge -- Percentage of city council seats elected at large

###################################
### Regression models
###################################
# only necessary if these packages are loaded
#detach("package:nlme")

## Look at simple correlations first
cor(data$mrp_ideology, data$policy, use="complete.obs")
cor(data$mrp_ideology, data$expenditures_capita, use="complete.obs")
cor(data$mrp_ideology, data$taxes_capita, use="complete.obs")
cor(data$mrp_ideology[data$local_salestax_allowed==1], data$salestax_share[data$local_salestax_allowed==1], use="complete.obs")


###################################
## Table 1:  Summary Statistics for Independent Variables
###################################

## Create table with summary statistics 
data$salestax_share[data$local_salestax_allowed!=1]<-NA
col1<-grep("mrp_ideology", colnames(data))
col2<-grep( "median_income",colnames(data))
col3<-grep( "city_pop",colnames(data))
col4<-grep( "percent_black",colnames(data))
col5<-grep( "house_value",colnames(data))
col6<-grep( "policy",colnames(data))
col7<-grep( "expenditures_capita",colnames(data))
col8<-grep( "taxes_capita",colnames(data))
col9<-grep( "salestax_share",colnames(data))

sum.data<-data[,c(col1, col2, col3, col4, col5, col6, col7, col8, col9)] 
sum.data<-sum.data[,-c(2,3,4)]
s<-t(summary(sum.data, digits=2) )
tab<-xtable(s, caption= "Summary Statistics of Key Variables", align=c("|c","|c","|c","|c","|c|","|c|","|c|","|c|")) 

###################################
## Table 2: Association Between City Liberalism and Policy Outcomes
###################################

##Baseline models with just city ideology.  
model0a<-(lmer(policy~mrp_ideology+(1|abb), data=data))
model0b<-(lmer(expenditures_capita~mrp_ideology+(1|abb), data=data))
model0c<-(lmer(taxes_capita~mrp_ideology+(1|abb), data=data))
model0d<-(lmer(salestax_share~mrp_ideology+(1|abb), data=data[data$local_salestax_allowed==1,]))

##Baseline models with city ideology and covariates for city demographics and income
model1a<-(lmer(policy~mrp_ideology+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
model1b<-(lmer(expenditures_capita~mrp_ideology+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
model1c<-(lmer(taxes_capita~mrp_ideology+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
model1d<-(lmer(salestax_share~mrp_ideology+median_income+city_pop+percent_black+house_value+(1|abb), data=data[data$local_salestax_allowed==1,]))

texreg(list(model0a,model1a, model0b,model1b,model0c,model1c,model0d,model1d), digits=2)

#stargazer(model0a,model1a, model0b,model1b,model0c,model1c,model0d,model1d, digits = 2, font.size="footnotesize",
#          title = "Association Between City Liberalism and Policy Outcomes",
#          covariate.labels = c("(Intercept)",   "Policy Conservatism",   "Median Income",  "City Population",  "Percent Black",   "Median Housing Value"),
#          dep.var.labels=c("Scaled Policy", "Per Capita Expend.","Per Capita Taxes","Sales Tax Share"),   omit.stat = "f", label = "models",column.sep.width="-1pt", star.cutoffs=c(.1,.05), intercept.top=T, intercept.bottom=F,digits.extra=0)



###################################
## Table 3: Effect of Institutions on Representation
###################################

model2a<-(lmer(policy~mrp_ideology*fog2+mrp_ideology*initiative2+mrp_ideology*partisan_elections2+mrp_ideology*term_limits2 +mrp_ideology*pal_binary+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
model2b<-(lmer(expenditures_capita~mrp_ideology*fog2+mrp_ideology*initiative2+mrp_ideology*partisan_elections2+mrp_ideology*term_limits2 +mrp_ideology*pal_binary+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
model2c<-(lmer(taxes_capita~mrp_ideology*fog2+mrp_ideology*initiative2+mrp_ideology*partisan_elections2+mrp_ideology*term_limits2 +mrp_ideology*pal_binary+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
model2d<-(lmer(salestax_share~mrp_ideology*fog2+mrp_ideology*initiative2+mrp_ideology*partisan_elections2+mrp_ideology*term_limits2 +mrp_ideology*pal_binary+median_income+city_pop+percent_black+house_value+(1|abb), data=data[data$local_salestax_allowed==1,]))

texreg(list(model0a,model2a, model0b,model2b,model0c,model2c,model0d,model2d))


#stargazer(model0a,model2a, model0b,model2b,model0c,model2c,model0d,model2d, digits = 2, font.size="footnotesize",
#          title = "Association Between City Liberalism and Policy Outcomes",
#          covariate.labels = c("(Intercept)", "Policy Conservatism","Elected Mayor","Direct Democracy","Partisan Elections" ,"Term Limits","At-large Elections"  ,"Median Income","City Population","Percent Black", "Median Housing Value","Conservatism x Mayor" ,"Conser. x Direct Dem.","Conser. x Part. Elect." ,"Conser. x Term Limits","Conser. x At-large"),
#          dep.var.labels=c("Scaled Policy", "Per Capita Expend.","Per Capita Taxes","Sales Tax Share"), omit.stat = "f", label = "models",column.sep.width="-1pt", star.cutoffs=c(.1,.05), intercept.top=T, intercept.bottom=F,digits.extra=0)


###################################
## Robustness check-  FE models (footnote 13)
###################################
model4a<-(lm(policy~mrp_ideology*fog2+mrp_ideology*initiative2+mrp_ideology*partisan_elections2+mrp_ideology*term_limits2 +mrp_ideology*pal_binary+median_income+city_pop+percent_black+house_value+factor(abb), data=data))
model4b<-(lm(expenditures_capita~mrp_ideology*fog2+mrp_ideology*initiative2+mrp_ideology*partisan_elections2+mrp_ideology*term_limits2 +mrp_ideology*pal_binary+median_income+city_pop+percent_black+house_value+factor(abb), data=data))
model4c<-(lm(taxes_capita~mrp_ideology*fog2+mrp_ideology*initiative2+mrp_ideology*partisan_elections2+mrp_ideology*term_limits2 +mrp_ideology*pal_binary+median_income+city_pop+percent_black+house_value+factor(abb), data=data))
model4d<-(lm(salestax_share~mrp_ideology*fog2+mrp_ideology*initiative2+mrp_ideology*partisan_elections2+mrp_ideology*term_limits2 +mrp_ideology*pal_binary+median_income+city_pop+percent_black+house_value+factor(abb), data=data[data$local_salestax_allowed==1,]))

texreg(list(model0a,model4a, model0b,model4b,model0c,model4c,model0d,model4d))

###################################
## Robustness check- RE models separately by institution (footnote 13)
###################################
#scaled policy 
summary(lmer(policy~mrp_ideology*fog2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(policy~mrp_ideology*initiative2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(policy~mrp_ideology*partisan_elections2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(policy~mrp_ideology*term_limits2+ median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(policy~mrp_ideology*pal_binary+ median_income+city_pop+percent_black+house_value+(1|abb), data=data))

#expenditures_capita
summary(lmer(expenditures_capita~mrp_ideology*fog2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(expenditures_capita~mrp_ideology*initiative2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(expenditures_capita~mrp_ideology*partisan_elections2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(expenditures_capita~mrp_ideology*term_limits2+ median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(expenditures_capita~mrp_ideology*pal_binary+ median_income+city_pop+percent_black+house_value+(1|abb), data=data))

#taxes_capita
summary(lmer(taxes_capita~mrp_ideology*fog2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(taxes_capita~mrp_ideology*initiative2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(taxes_capita~mrp_ideology*partisan_elections2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(taxes_capita~mrp_ideology*term_limits2+ median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(taxes_capita~mrp_ideology*pal_binary+ median_income+city_pop+percent_black+house_value+(1|abb), data=data))

#salestax_share
summary(lmer(salestax_share~mrp_ideology*fog2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(salestax_share~mrp_ideology*initiative2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(salestax_share~mrp_ideology*partisan_elections2+median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(salestax_share~mrp_ideology*term_limits2+ median_income+city_pop+percent_black+house_value+(1|abb), data=data))
summary(lmer(salestax_share~mrp_ideology*pal_binary+ median_income+city_pop+percent_black+house_value+(1|abb), data=data))



##########################
## VALIDATION AND DESCRIPTIVE GRAPHS
#########################

## correlation between ideology and presidential vote share in the paper
cor(data$mrp_ideology, data$pres_2008, use="complete.obs")

###########
## Figure 1: Policy preferences of large cities
#############
pdf(paste("figure1_cities_large.pdf",sep=""),width=7.5, height=10)
data2<-subset(data, city_pop>2.5)
 par(mfrow=c(1,1),mar=c(3,10,3,3))
data.order<-order(as.numeric(data2$mrp_ideology), decreasing = FALSE)
data2<-data2[data.order,]
right.label<-paste("(", round(data2$mrp_ideology,2), ")", sep="")
left.label<-data2$city
left.label<-paste("",left.label,", ",as.vector(data2$abb), sep="")
left.label<-gsub(" (balance) ", "", left.label,fixed = T)
left.label<-gsub("(balance) ", "", left.label,fixed = T)
left.label<-gsub("(balance)", "", left.label,fixed = T)
left.label<-gsub(" ,", ",", left.label,fixed = T)
y.axis<-c(1: dim(data2)[1])
plot(data2$mrp_ideology, y.axis, axes=F, xlim=c(-1.2,1),ylim=c(.1,dim(data2)[1]),pch=19,main="",xlab='',ylab='',type="n")
for (i in 1: dim(data2)[1]){
	 segments(-3, i, 3, i, col="grey")
}
segments(0,0,0,(dim(data2)[1]+1), col="light grey")
segments(data2$mrp_ideology_lower, y.axis+.2, data2$mrp_ideology_upper, y.axis+.2, col="black",lwd=2)
points(data2$mrp_ideology, y.axis+.2, col="black", pch=19, cex=.75)
axis(1, at = c(-3,-2.5,-2,-1.5,-1,-.5,0,.5,1,1.5,2,2.5,3), pos=0,cex.axis=.7 )
axis(2, at = c(0,y.axis, (dim(data2)[1]+2)), label = c('',left.label,''), las = 1, cex.axis = 1, mgp =
c(2,1,0), pos=-1.2,cex.axis=.75, tick=FALSE)
axis(4,at = c(0,y.axis, (dim(data2)[1]+2)), pos=1,label = c('',right.label,''), tick=FALSE,las = 2,cex.axis=.75)
title(main="Policy Preferences of Mass Public by City
(More than 250,000 people)",xlab = "", ylab = "", cex.lab=1)
title(xlab = "City Conservatism",cex.lab=1, line=1)
dev.off()

###########
## Figures 2 and 3: Bivariate analysis
#############

dependentvar <- c("policy","salestax_share","taxes_capita","expenditures_capita")
depvarnames <- c("Policy Scale","Share of Taxes From Sales Taxes","Taxes per Capita","Expenditures per Capita")

bivariateplot <- function(dependentvar,depvarnames,abl=F,span=.75,...){
    attach(data)
    for (i in dependentvar){
       dvar <- get(i)
       colors <- ifelse(term_limits2==1,"red","black")
       colors[is.na(colors)] <- "blue"
      par(mar = c(5,6,4,2) + 0.1)
plot(dvar~mrp_ideology,axes=F,ylab="",xlab="City Policy Conservatism",cex=log(city_pop/mean(city_pop,na.rm=T)*2),col="grey",...)
  title(ylab =depvarnames[dependentvar==i], cex.lab = 1,
      line = 3.25)
    axis(1)
       axis(2, las=2)
       temp <- data.frame(x=mrp_ideology, y=dvar)
       y.loess <- loess(y ~ x,temp,span=span)
       y.predict <- predict(y.loess, data.frame(x=seq(min(temp$x,na.rm=T),max(temp$x,na.rm=T),sd(temp$x,na.rm=T)/5)))
       lines(y.predict~seq(min(temp$x,na.rm=T),max(temp$x,na.rm=T),sd(temp$x,na.rm=T)/5),col="black")
       if(abl==T){abline(lm(temp$y~temp$x),col="black")}
       
       
       ind <- which(city=="Los Angeles")
       text(mrp_ideology[ind],dvar[ind],"LA",col="black")
       ind <- which(city=="New York")
       text(mrp_ideology[ind],dvar[ind],"NY",col="black")
       ind <- which(city=="Chicago")
       text(mrp_ideology[ind],dvar[ind],"CHI",col="black")
       ind <- which(city=="Phoenix")
       text(mrp_ideology[ind],dvar[ind],"PHX",col="black")
       ind <- which(city=="Detroit")
       text(mrp_ideology[ind],dvar[ind],"DET",col="black")
       ind <- which(city=="San Diego")
       text(mrp_ideology[ind],dvar[ind],"SD",col="black")
       ind <- which(city=="Philadelphia")
       text(mrp_ideology[ind],dvar[ind],"PHI",col="black")
       ind <- which(city=="Dallas")
       text(mrp_ideology[ind],dvar[ind],"DAL",col="black")
       ind <- which(city=="Houston")
       text(mrp_ideology[ind],dvar[ind],"HOU",col="black")
       ind <- which(city=="San Antonio")
       text(mrp_ideology[ind],dvar[ind],"SA",col="black")
       ind <- which(city=="Washington")
       text(mrp_ideology[ind],dvar[ind],"DC",col="black")
    }
  }

#############
# Figure 2
#############
pdf("figure2_responsiveness.pdf",height=10,width=10)
par(mfcol=c(2,2))
bivariateplot(dependentvar,depvarnames,span=1)
dev.off()

#############
# Figure 3
#############
pdf("figure3_zoom.pdf",height=5,width=10)
par(mfcol=c(1,2))
dependentvarTemp <- c("taxes_capita")
depvarnamesTemp <- c("Taxes per Capita")
bivariateplot(dependentvarTemp,depvarnamesTemp,abl=T,ylim=c(0,2200))
dependentvarTemp <- c("expenditures_capita")
depvarnamesTemp <- c("Expenditures per Capita")
bivariateplot(dependentvarTemp,depvarnamesTemp,abl=T,ylim=c(0,3200))
dev.off()


###########
## Figure 4: Bayesian model error distribution
## Note: This section is computationally intensive
#############

dependentvar <- c("policy","salestax_share","taxes_capita","expenditures_capita")
depvarnames <- c("Policy Scale","Share of Taxes From Sales Taxes","Taxes per Capita","Expenditures per Capita")


data2<-subset(data, !is.na(house_value))
modelfile <- "regression.jags"

bayesResidual <- function(dependentvar,depvarnames){
  MAElist <- list()
  attach(data2)
    
  for (i in dependentvar){
     iterations=1000
     thin=1
     adapt=100

     jdata <- list(y=get(i),
                   x=cbind(mrp_ideology,
                           rep(NA,dim(data2)[1]),
                           rep(NA,dim(data2)[1]),
                           rep(NA,dim(data2)[1]),
                           rep(NA,dim(data2)[1]),
                           rep(NA,dim(data2)[1]),
                           median_income,
                           city_pop,
                           initiative2,
                           partisan_elections2,
                           term_limits2,
                           fog2,
                           pal_binary,
                           percent_black,
                           house_value),
                   k=15,n=dim(data2)[1])

     print("Preparing the model...")
     system.time(
     model <- jags.model(file=modelfile,data=jdata,
                         n.adapt=adapt))
     print("Running the model...")
     system.time(
     output <- coda.samples(model,
                            variable.names=c("alpha","beta","x"),
                            n.iter=iterations,thin=thin))
     coefnames <- c("(Intercept)","mrp_ideology","mrp_ideology:initiative2","mrp_ideology:partisan_elections2","mrp_ideology:term_limits2","mrp_ideology:fog2","mrp_ideology:pal_binary","median_income","city_pop","initiative2","partisan_elections2","term_limits2","fog2","pal_binary","percent_black","house_value")

     estimatedxs <- matrix(colMeans(output[[1]][,grep("x",colnames(output[[1]]))]),jdata$n,jdata$k)

     colnames(output[[1]][,1:length(coefnames)]) <- coefnames
     outlist <- list(output=output[[1]],jdata=jdata,xs=estimatedxs)

     ymatrix <- matrix(rep(outlist$jdata$y,dim(output[[1]])[1]),dim(output[[1]])[1],outlist$jdata$n,byrow=T)[,-which(is.na(get(i)))]
     MAE <- abs(ymatrix - output[[1]][,1:length(coefnames)] %*% t(cbind(1,estimatedxs[-which(is.na(get(i))),])))

     hist(MAE[MAE < quantile(MAE,.975)],axes=F,main=paste("",depvarnames[dependentvar==i],sep=""),ylab="",xlab="Error")
     axis(1)

     MAElist[[which(dependentvar==i)]] <- MAE
   }
  return(MAElist)
}

#############
# Figure 4 
#############
pdf(file="figure4_ErrorDistribution2.pdf")
par(mfcol=c(2,2))
MAElist <- bayesResidual(dependentvar,depvarnames)
dev.off()

# median error numbers used in the text
median(MAElist[[1]])
median(MAElist[[2]])
median(MAElist[[3]])
median(MAElist[[4]])

###########
## Figures 5 through 9: Estimating Institutional
##     Effects with Entropy Balancing
#############

data2<-subset(data, !is.na(data$partisan_elections2))
data2<-subset(data2, !is.na(data2$initiative2))
data2<-subset(data2, !is.na(data2$median_income))
data2<-subset(data2, !is.na(data2$city_pop))
data2<-subset(data2, !is.na(data2$house_value))
data2<-subset(data2, !is.na(data2$fog2)) 
data2<-subset(data2, !is.na(data2$term_limits2))
data2<-subset(data2, !is.na(data2$pal_binary))
data2<-subset(data2, !is.na(data2$mrp_ideology))

dependentvar <- c("policy","expenditures_capita","taxes_capita","salestax_share")
depvarnames <- c("Policy Scale","Expenditures per Capita","Taxes per Capita","Share of Taxes From Sales Taxes")

draw_loess_ci <- function(y,x,w,span,iter=1000,...){
       y.loess <- loess(y ~ x,span=span,weights=w)

       fakex <- seq(min(x,na.rm=T),max(x,na.rm=T),sd(x,na.rm=T)/5)
       y.predict <- predict(y.loess, data.frame(x=fakex))
       
       lines(y.predict~fakex,lwd=3, ...)
} 

matchingplot <- function(dependentvar,depvarnames,datname,institution,span=.75,...){
  j<-1
    for (i in dependentvar){
       dat <- get(datname)
       print(i)
       dvar <- dat[,i]
       dat<-subset(get(datname), ! is.na(dvar))
       if(i=="salestax_share"){dat <- dat[dat$local_salestax_allowed==1,]}
       dvar <- dat[,i]
       inst <- dat[[institution]]
       inst <- dat[[institution]]
       colors <- ifelse(inst==1,"black","grey")
       colors[is.na(colors)] <- "white"

       # do entropy balancing
       institutions <- c("fog2","initiative2","partisan_elections2","term_limits2","pal_binary")
       covars<-cbind(dat[,names(dat) %in% institutions[institutions != institution]], dat$median_income, dat$city_pop, dat$house_value, dat$mrp_ideology,  dat$mrp_ideology^2)
       ebweights <- ebalance(Treatment=dat[[institution]],X=covars)$w
       dat$w <- 1
       dat$w[dat[[institution]]==0] <- ebweights

       ### DRAW THE GRAPH
   par(mar = c(5,6,4,2) + 0.1) 
    plot(dvar~dat[["mrp_ideology"]],axes=F,xlim=c(-1,.5),
       ylab="",xlab="City Policy Conservatism",col=colors,pch=1,...)
  title(ylab =depvarnames[dependentvar==i], cex.lab = 1,
      line = 3.25)
      
       axis(1,at=c(-1,-.5,0,.5))
       axis(2,las=2)

       temp <- data.frame(x=dat[["mrp_ideology"]][inst==1], y=dvar[inst==1],w=dat$w[inst==1])
       draw_loess_ci(temp$y,temp$x,temp$w,span,col="black")
       temp <- data.frame(x=dat[["mrp_ideology"]][inst==0], y=dvar[inst==0],w=dat$w[inst==0])
       draw_loess_ci(temp$y,temp$x,temp$w,span,col="grey")
		dat$mrp_ideology.sq<-dat$mrp_ideology*dat$mrp_ideology

       ### PUT TOGETHER THE REGRESSION TABLE
      if (j==1){
       mat1<-lm(dvar~dat[["mrp_ideology"]]*inst+median_income+city_pop+house_value+factor(abb),weights=w, data=dat)
         	assign("mat1",mat1,envir = .GlobalEnv)
  	}
        if (j==2){
       mat2<-lm(dvar~dat[["mrp_ideology"]]*inst+median_income+city_pop+house_value+factor(abb),weights=w, data=dat)
          	assign("mat2",mat2,envir = .GlobalEnv)
 	}      
    	if (j==3){
       mat3<-lm(dvar~dat[["mrp_ideology"]]*inst+median_income+city_pop+house_value+factor(abb),weights=w, data=dat)
       	assign("mat3",mat3,envir = .GlobalEnv)
    	}      
    	if (j==4){
       mat4<-lm(dvar~dat[["mrp_ideology"]]*inst+median_income+city_pop+house_value+factor(abb),weights=w, data=dat)
       	assign("mat4",mat4,envir = .GlobalEnv)
    	}  	
    j<- j+1

 }
 print(texreg(list(mat1, mat2,  mat3, mat4)))
}


#############
# Figure 5
#############
pdf("figure5_EBALANCEgov2.pdf")
par(mfrow=c(2,2))
institution <- "fog2"
matchingplot(dependentvar,depvarnames,datname="data2",institution,span=.85,cex=.5)
dev.off()

#############
# Figure 6
#############
pdf("figure6_EBALANCEpartisan_elections2.pdf")
par(mfrow=c(2,2))
institution <- "partisan_elections2"
matchingplot(dependentvar,depvarnames,datname="data2",institution,span=.85)
dev.off()

#############
# Figure 7
#############
pdf("figure7_EBALANCEinitiative2.pdf")
par(mfrow=c(2,2))
institution <- "initiative2"
matchingplot(dependentvar,depvarnames,datname="data2",institution,span=1,cex=.75)
dev.off()

#############
# Figure 8
#############
pdf("figure8_EBALANCEterm_limits2.pdf")
par(mfrow=c(2,2))
institution <- "term_limits2"
matchingplot(dependentvar,depvarnames,datname="data2",institution,span=1,cex=.5)
dev.off()

#############
# Figure 9
#############
pdf("figure9_EBALANCEpal_binary.pdf")
par(mfrow=c(2,2))
institution <- "pal_binary"
matchingplot(dependentvar,depvarnames,datname="data2",institution,span=1,cex=.5)
dev.off()


###################
## Figure 10: Policy preferences in 4 states
###################
library(mgcv)

pdf(paste("figure10_cities_states.pdf",sep=""),width=7.5, height=10)
 par(mfrow=c(1,1),mar=c(3,6.5,3,2))
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
data2<-subset(data, city_pop>.75 & abb=="TX")
data.order<-order(as.numeric(data2$mrp_ideology), decreasing = FALSE)
data2<-data2[data.order,]
right.label<-paste("(", round(data2$mrp_ideology,2), ")", sep="")
left.label<-data2$city
left.label<-paste("",left.label,", ",as.vector(data2$abb), sep="")
left.label<-gsub(" (balance) ", "", left.label,fixed = T)
left.label<-gsub("(balance) ", "", left.label,fixed = T)
left.label<-gsub("(balance)", "", left.label,fixed = T)
left.label<-gsub(" ,", ",", left.label,fixed = T)
y.axis<-c(1: dim(data2)[1])
plot(data2$mrp_ideology, y.axis, axes=F, xlim=c(-1.2,1),ylim=c(.1,dim(data2)[1]),pch=19,main="",xlab='',ylab='',type="n")

for (i in 1: dim(data2)[1]){
	 segments(-3, i, 3, i, col="grey")
}
segments(0,0,0,(dim(data2)[1]+1), col="light grey")
segments(data2$mrp_ideology_lower, y.axis+.2, data2$mrp_ideology_upper, y.axis+.2, col="black",lwd=2)
points(data2$mrp_ideology, y.axis+.2, col="black", pch=19, cex=.75)
axis(1, at = c(-3,-2.5,-2,-1.5,-1,-.5,0,.5,1,1.5,2,2.5,3), pos=0,cex.axis=.7 )
axis(2, at = c(0,y.axis, (dim(data2)[1]+2)), label = c('',left.label,''), las = 1, cex.axis = 1, mgp =
c(2,1,0), pos=-1.2,cex.axis=.75, tick=FALSE)
title(main="City Policy Preferences in Texas",  cex.main=.9, line=1)
title(xlab = "City Conservatism",cex.lab=1, line=1)

data2<-subset(data, city_pop>.75 & abb=="VA")
data.order<-order(as.numeric(data2$mrp_ideology), decreasing = FALSE)
data2<-data2[data.order,]
right.label<-paste("(", round(data2$mrp_ideology,2), ")", sep="")
left.label<-data2$city
left.label<-paste("",left.label,", ",as.vector(data2$abb), sep="")
left.label<-gsub(" (balance) ", "", left.label,fixed = T)
left.label<-gsub("(balance) ", "", left.label,fixed = T)
left.label<-gsub("(balance)", "", left.label,fixed = T)
left.label<-gsub(" ,", ",", left.label,fixed = T)
y.axis<-c(1: dim(data2)[1])
plot(data2$mrp_ideology, y.axis, axes=F, xlim=c(-1.2,1),ylim=c(.1,dim(data2)[1]),pch=19,main="",xlab='',ylab='',type="n")
for (i in 1: dim(data2)[1]){
	 segments(-3, i, 3, i, col="grey")
}
segments(0,0,0,(dim(data2)[1]+1), col="light grey")
segments(data2$mrp_ideology_lower, y.axis+.2, data2$mrp_ideology_upper, y.axis+.2, col="black",lwd=2)
points(data2$mrp_ideology, y.axis+.2, col="black", pch=19, cex=.75)
axis(1, at = c(-3,-2.5,-2,-1.5,-1,-.5,0,.5,1,1.5,2,2.5,3), pos=0,cex.axis=.7 )
axis(2, at = c(0,y.axis, (dim(data2)[1]+2)), label = c('',left.label,''), las = 1, cex.axis = 1, mgp =
c(2,1,0), pos=-1.2,cex.axis=.75, tick=FALSE)
title(main="City Policy Preferences in Virginia",xlab = "", ylab = "", cex.main=.9)
title(xlab = "City Conservatism",cex.lab=1, line=1)

data2<-subset(data, city_pop>.75 & abb=="MI")
data.order<-order(as.numeric(data2$mrp_ideology), decreasing = FALSE)
data2<-data2[data.order,]
right.label<-paste("(", round(data2$mrp_ideology,2), ")", sep="")
left.label<-data2$city
left.label<-paste("",left.label,", ",as.vector(data2$abb), sep="")
left.label<-gsub(" (balance) ", "", left.label,fixed = T)
left.label<-gsub("(balance) ", "", left.label,fixed = T)
left.label<-gsub("(balance)", "", left.label,fixed = T)
left.label<-gsub(" ,", ",", left.label,fixed = T)
y.axis<-c(1: dim(data2)[1])
plot(data2$mrp_ideology, y.axis, axes=F, xlim=c(-1.2,1),ylim=c(.1,dim(data2)[1]),pch=19,main="",xlab='',ylab='',type="n")
for (i in 1: dim(data2)[1]){
	 segments(-3, i, 3, i, col="grey")
}
segments(0,0,0,(dim(data2)[1]+1), col="light grey")
segments(data2$mrp_ideology_lower, y.axis+.2, data2$mrp_ideology_upper, y.axis+.2, col="black",lwd=2)
points(data2$mrp_ideology, y.axis+.2, col="black", pch=19, cex=.75)
axis(1, at = c(-3,-2.5,-2,-1.5,-1,-.5,0,.5,1,1.5,2,2.5,3), pos=0,cex.axis=.7 )
axis(2, at = c(0,y.axis, (dim(data2)[1]+2)), label = c('',left.label,''), las = 1, cex.axis = 1, mgp =
c(2,1,0), pos=-1.2,cex.axis=.75, tick=FALSE)
title(main="City Policy Preferences in Michigan",xlab = "", ylab = "", cex.main=.9)
title(xlab = "City Conservatism",cex.lab=1, line=1)

data2<-subset(data, city_pop>.75 & abb=="MA")
data.order<-order(as.numeric(data2$mrp_ideology), decreasing = FALSE)
data2<-data2[data.order,]
right.label<-paste("(", round(data2$mrp_ideology,2), ")", sep="")
left.label<-data2$city
left.label<-paste("",left.label,", ",as.vector(data2$abb), sep="")
left.label<-gsub(" (balance) ", "", left.label,fixed = T)
left.label<-gsub("(balance) ", "", left.label,fixed = T)
left.label<-gsub("(balance)", "", left.label,fixed = T)
left.label<-gsub(" ,", ",", left.label,fixed = T)
y.axis<-c(1: dim(data2)[1])
plot(data2$mrp_ideology, y.axis, axes=F, xlim=c(-1.2,1),ylim=c(.1,dim(data2)[1]),pch=19,main="",xlab='',ylab='',type="n")
for (i in 1: dim(data2)[1]){
	 segments(-3, i, 3, i, col="grey")
}
segments(0,0,0,(dim(data2)[1]+1), col="light grey")
segments(data2$mrp_ideology_lower, y.axis+.2, data2$mrp_ideology_upper, y.axis+.2, col="black",lwd=2)
points(data2$mrp_ideology, y.axis+.2, col="black", pch=19, cex=.75)
axis(1, at = c(-3,-2.5,-2,-1.5,-1,-.5,0,.5,1,1.5,2,2.5,3), pos=0,cex.axis=.7 )
axis(2, at = c(0,y.axis, (dim(data2)[1]+2)), label = c('',left.label,''), las = 1, cex.axis = 1, mgp =
c(2,1,0), pos=-1.2,cex.axis=.75, tick=FALSE)
title(main="City Policy Preferences in Massachusetts",xlab = "", ylab = "", cex.main=.9)
title(xlab = "City Conservatism",cex.lab=1, line=1)
dev.off()

# Get State names
states<-read.csv(file="states.csv")
states$state<-gsub(" ","",states$state)
states$state[states$state=="DistrictofColumbia"]<-"District of Columbia"
states$state[states$state=="NewHampshire"]<-"New Hampshire"
states$state[states$state=="NewJersey"]<-"New Jersey"
states$state[states$state=="NewMexico"]<-"New Mexico"
states$state[states$state=="NewYork"]<-"New York"
states$state[states$state=="NorthCarolina"]<-"North Carolina"
states$state[states$state=="NorthDakota"]<-"North Dakota"
states$state[states$state=="RhodeIsland"]<-"Rhode Island"
states$state[states$state=="SouthCarolina"]<-"South Carolina"
states$state[states$state=="SouthDakota"]<-"South Dakota"
states$state[states$state=="WestVirginia"]<-"West Virginia"

data$city_pop<-data$city_pop*100000
knight_survey_2002<-read.csv("knight_survey_2002.csv")
data<-merge(data, knight_survey_2002, by="city_id", all.x=T)
bayarea_conservatism<-read.csv("cities_lib.csv")
data<-merge(data, states, by.x=c("abb" ), by.y=c("abb"), all.x=T)
data<-merge(data, bayarea_conservatism, by.x=c("state", "city.x"), by.y=c("state", "city"), all.x=T)
cor(data$mrp_ideology, data$conservatism, use="complete.obs")
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

data$pres_2008_2<-as.numeric(as.vector(data$pres_2008))
data$city_pop<-as.numeric(as.vector(data$city_pop))
data$pres_2008_2<-1-data$pres_2008_2
data$self_ideology_knight_2<-7-data$self_ideology_knight

cors<-rep(NA, 4)
cors[1]<-cor(data$mrp_ideology, data$raw_ideology, use="complete.obs")
cors[2]<-cor(data$mrp_ideology[data$city_pop > 100000], data$raw_ideology[data$city_pop > 100000], use="complete.obs")
cors[3]<-cor(data$mrp_ideology, data$pres_2008_2, use="complete.obs")
cors[4]<-cor(data$mrp_ideology, data$self_ideology_knight_2, use="complete.obs")
cor(data$mrp_ideology, data$conservatism, use="complete.obs")
data$conservatism<-data$conservatism/100
data$city_pop<-data$city_pop/100000

p1<-(ggplot(data=data, aes(x=mrp_ideology, y=raw_ideology))
          + geom_point()
          + opts(axis.title.x=theme_text(size=12))
		  + opts(axis.title.y=theme_text(size=12))
			+stat_smooth(se=TRUE,method = "bam", formula =y ~ s(x, k = 9))
        + ylab("Raw Policy Conservatism (all cities)")
        + xlab("MRP-based Estimate of City Conservatism"))	

p2<-(ggplot(data=data[data$city_pop > 1,], aes(x=mrp_ideology, y=raw_ideology))
          + geom_point()
          + opts(axis.title.x=theme_text(size=12))
		  + opts(axis.title.y=theme_text(size=12))
			+stat_smooth(se=TRUE,method = "bam", formula =y ~ s(x, k = 9))
        + ylab("Raw Policy Conservatism (large cities)")
        + xlab("MRP-based Estimate of City Conservatism"))	

p3<-(ggplot(data=data, aes(x=mrp_ideology, y=pres_2008_2))
          + geom_point()
          + opts(axis.title.x=theme_text(size=12))
		  + opts(axis.title.y=theme_text(size=12))
			+stat_smooth(se=TRUE,method = "bam", formula =y ~ s(x, k = 9))
        + ylab("Republican Presidential Vote (2008)")
        + xlab("MRP-based Estimate of City Conservatism"))	

p4<-(ggplot(data=data, aes(x=mrp_ideology, y=self_ideology_knight_2))
          + geom_point()
          + opts(axis.title.x=theme_text(size=12))
		  + opts(axis.title.y=theme_text(size=12))
		+stat_smooth(se=TRUE,method = "bam", formula =y ~ s(x, k = 9))
        + ylab("Conservatism of 24 Cities from Knight Survey")
        + xlab("MRP-based Estimate of City Conservatism"))	
     
p5<-(ggplot(data=data, aes(x=mrp_ideology, y=conservatism))
          + geom_point()
          + opts(axis.title.x=theme_text(size=12))
		  + opts(axis.title.y=theme_text(size=12))
		+stat_smooth(se=TRUE,method = "bam", formula =y ~ s(x, k = 9))
        + ylab("Republican Presidential Vote (2004)")
        + xlab("MRP-based Estimate of City Conservatism"))	
        
               
#####
### Figure 11: Relationship between MRP and Disaggregated City Policy Conservatism
####
pdf(file="figure11_validation_internal.pdf", width=7, height=8)
multiplot(p1, p2,  cols=1)
dev.off()

#####
### Figure 12: Relationship between City Policy Conservatism and Presidential Vote Share
####
pdf(file="figure12_validation_pvote.pdf", width=7, height=8)
multiplot(p3, p5, cols=1)
dev.off()
    
#####
### Figure 13: Relationship between City Policy Conservatism and External Metric of City Conservatism
####
pdf(file="figure13_validation_external.pdf", width=7, height=8)
multiplot(p4, cols=1)
dev.off()
    
########    
### Figure 14: Dimensionality
########

load("cityFederal.RData")

outputFederal$xbar[,1] <- -outputFederal$xbar[,1]

cor(outputCity$xbar[,1],outputFederal$xbar[,1])

pdf("figure14_municipalFederalCorrelation.pdf")
plot(outputCity$xbar[,1]~outputFederal$xbar[,1],axes=F,main="",ylab="Municipal policy scale",xlab="Federal policy scale")
axis(1)
axis(2)
abline(lm(outputCity$xbar[,1]~outputFederal$xbar[,1]))
text(-2,1.5,"Correlation=0.75")
dev.off()

#############
## CEM Matching model
############
library(cem)

### remove extraneous variables
data<-data[,-which(colnames(data)== "perc_atlarge2")]
data<-data[,-which(colnames(data)== "raw_ideology")]
data<-data[,-which(colnames(data)== "self_ideology_knight")]
data<-data[,-which(colnames(data)== "sample_size")]
data<-data[,-which(colnames(data)== "county")]
data<-data[,-which(colnames(data)== "region")]
data<-data[,-which(colnames(data)== "conservatism")]
data<-data[,-which(colnames(data)== "self_ideology_knight_2")]
data<-data[,-which(colnames(data)== "pres_2008_2")]
data<-data[,-which(colnames(data)== "city.y")]
data<-data[,-which(colnames(data)== "pop")]
data<-data[,-which(colnames(data)== "region_political")]
data<-data[,-which(colnames(data)== "mrp_ideology_sd")]
data<-data[,-which(colnames(data)== "mrp_ideology_lower")]
data<-data[,-which(colnames(data)== "mrp_ideology_upper")]


institutions<-  list(1, 0)

quantile(data$house_value, probs = seq(0, 1, 0.2),  na.rm = T)
quantile(data$city_pop, probs = seq(0, 1, 0.2),  na.rm = T)
quantile(data$median_income, probs = seq(0, 1, 0.2),  na.rm = T)
quantile(data$percent_black, probs = seq(0, 1, 0.2),  na.rm = T)
quantile(data$mrp_ideology, probs = seq(0, 1, 0.2),  na.rm = T)

## Tell CEM the cutpoints to form the coursened exact matches for each of our key variables.  In this case, I used quantiles (0, 20, 40, 60, 80, 100).  But there are probably other equally sensible approaches.
house_cut <- c(0, quantile(data$house_value, probs = seq(0, 1, 0.2),  na.rm = T)[2],quantile(data$house_value, probs = seq(0, 1, 0.2),  na.rm = T)[3],quantile(data$house_value, probs = seq(0, 1, 0.2),  na.rm = T)[4],quantile(data$house_value, probs = seq(0, 1, 0.2),  na.rm = T)[5], 1000)
pop_cut <- c(0, quantile(data$city_pop, probs = seq(0, 1, 0.2),  na.rm = T)[2],quantile(data$city_pop, probs = seq(0, 1, 0.2),  na.rm = T)[3],quantile(data$city_pop, probs = seq(0, 1, 0.2),  na.rm = T)[4],quantile(data$city_pop, probs = seq(0, 1, 0.2),  na.rm = T)[5],3,10, 100)
income_cut <- c(0, quantile(data$median_income, probs = seq(0, 1, 0.2),  na.rm = T)[2],quantile(data$median_income, probs = seq(0, 1, 0.2),  na.rm = T)[3],quantile(data$median_income, probs = seq(0, 1, 0.2),  na.rm = T)[4],quantile(data$median_income, probs = seq(0, 1, 0.2),  na.rm = T)[5], 1000)
black_cut<-c(0, quantile(data$percent_black, probs = seq(0, 1, 0.2),  na.rm = T)[2],quantile(data$percent_black, probs = seq(0, 1, 0.2),  na.rm = T)[3],quantile(data$percent_black, probs = seq(0, 1, 0.2),  na.rm = T)[4],quantile(data$percent_black, probs = seq(0, 1, 0.2),  na.rm = T)[5], 1)
mrp_cut<-c(-2, quantile(data$mrp_ideology, probs = seq(0, 1, 0.2),  na.rm = T)[2],quantile(data$mrp_ideology, probs = seq(0, 1, 0.2),  na.rm = T)[3],quantile(data$mrp_ideology, probs = seq(0, 1, 0.2),  na.rm = T)[4],quantile(data$mrp_ideology, probs = seq(0, 1, 0.2),  na.rm = T)[5], 2)

## Tell CEM the cutpoints to form the coursened exact matches for each of our key variables.  In this case, I used quantiles (0, 25, 50, 75,100).  But there are probably other equally sensible approaches.
house_cut2 <- c(0, quantile(data$house_value, probs = seq(0, 1, 0.25),  na.rm = T)[2],quantile(data$house_value, probs = seq(0, 1, 0.25),  na.rm = T)[3],quantile(data$house_value, probs = seq(0, 1, 0.25),  na.rm = T)[4], 1000)
pop_cut2 <- c(0, quantile(data$city_pop, probs = seq(0, 1, 0.25),  na.rm = T)[2],quantile(data$city_pop, probs = seq(0, 1, 0.25),  na.rm = T)[3],quantile(data$city_pop, probs = seq(0, 1, 0.25),  na.rm = T)[4],3,10, 100)
income_cut2 <- c(0, quantile(data$median_income, probs = seq(0, 1, 0.25),  na.rm = T)[2],quantile(data$median_income, probs = seq(0, 1, 0.25),  na.rm = T)[3],quantile(data$median_income, probs = seq(0, 1, 0.25),  na.rm = T)[4], 1000)
black_cut2<-c(0, quantile(data$percent_black, probs = seq(0, 1, 0.25),  na.rm = T)[2],quantile(data$percent_black, probs = seq(0, 1, 0.25),  na.rm = T)[3],quantile(data$percent_black, probs = seq(0, 1, 0.25),  na.rm = T)[4],1)
mrp_cut2<-c(-2, quantile(data$mrp_ideology, probs = seq(0, 1, 0.25),  na.rm = T)[2],quantile(data$mrp_ideology, probs = seq(0, 1, 0.25),  na.rm = T)[3],quantile(data$mrp_ideology, probs = seq(0, 1, 0.25),  na.rm = T)[4], 2)


## Run four sets of models. One for each institution.


#####
## gov type (#1)
#####

##Policy
data2<-subset(data, !is.na(policy))
mat2a <- cem(treatment = "fog2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita",  "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",   "fog2",  "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat2a$breaks
data2_model2<-data2
data2_model2$match<-mat2a$matched
data2_model2$weight<-mat2a$w
data2_model2<-subset(data2_model2, match=="TRUE")

summary(att(mat2a,policy~mrp_ideology*fog2 +median_income+percent_black+city_pop+house_value+factor(abb), data=data2))
fog_policy<-lm(policy~mrp_ideology*fog2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model2)

##expenditures_capita
data2<-subset(data, !is.na(expenditures_capita))
mat2b <- cem(treatment = "fog2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita",  "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",   "fog2", "fog2", "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat2b$breaks
data2_model2<-data2
data2_model2$match<-mat2b$matched
data2_model2$weight<-mat2b$w
data2_model2<-subset(data2_model2, match=="TRUE")

summary(att(mat2b,expenditures_capita~mrp_ideology*fog2+median_income+percent_black+city_pop+house_value+factor(abb), data=data2))
fog_exp_pc<-lm(expenditures_capita~mrp_ideology*fog2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model2)

##taxes_capita
data2<-subset(data, !is.na(taxes_capita))
mat2c <- cem(treatment = "fog2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita",  "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",   "fog2", "fog2", "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat2c$breaks
data2_model2<-data2
data2_model2$match<-mat2c$matched
data2_model2$weight<-mat2c$w
data2_model2<-subset(data2_model2, match=="TRUE")

summary(att(mat2c,taxes_capita~mrp_ideology*fog2+median_income+percent_black+city_pop+house_value+factor(abb), data=data2))
fog_taxes_pc<-lm(taxes_capita~mrp_ideology*fog2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model2)

##salestax_share
data2<-subset(data, !is.na(salestax_share) & local_salestax_allowed==1 )
mat2d.st <- cem(treatment = "fog2", data = data2[data2$local_salestax_allowed==1,], drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita",  "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita", "fog2",  "fog2",  "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat2c$breaks
data2_model2<-data2
data2_model2$match<-mat2d.st$matched
data2_model2$weight<-mat2d.st$w
data2_model2<-subset(data2_model2, match=="TRUE")

summary(att(mat2d.st,salestax_share~mrp_ideology*partisan_elections2 +median_income+percent_black+city_pop+house_value+factor(abb), data=data2[data2$local_salestax_allowed==1,]))

fog_salestax<-lm(salestax_share~mrp_ideology*fog2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model2[data2_model2$local_salestax_allowed==1,])


#####
## partisan_elections (#2)
#####

##Policy
data2<-subset(data, !is.na(policy))
mat3a <- cem(treatment = "partisan_elections2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita",  "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "partisan_elections2",   "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat3a$breaks
data2_model3<-data2
data2_model3$match<-mat3a$matched
data2_model3$weight<-mat3a$w
data2_model3<-subset(data2_model3, match=="TRUE")

summary(att(mat3a,policy~mrp_ideology*partisan_elections2+percent_black +median_income+percent_black+city_pop+house_value+factor(abb), data=data2))
partisanelections_policy<-lm(policy~mrp_ideology*partisan_elections2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model3)

##expenditures_capita
data2<-subset(data, !is.na(expenditures_capita))
mat3b <- cem(treatment = "partisan_elections2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita",  "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "partisan_elections2",  "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat3b$breaks
data2_model3<-data2
data2_model3$match<-mat3b$matched
data2_model3$weight<-mat3b$w
data2_model3<-subset(data2_model3, match=="TRUE")

summary(att(mat3b,expenditures_capita~mrp_ideology*partisan_elections2+median_income+percent_black+city_pop+house_value+factor(abb), data=data2))
partisanelections_exp_pc<-lm(expenditures_capita~mrp_ideology*partisan_elections2+median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model3)

##taxes_capita
data2<-subset(data, !is.na(taxes_capita))
mat3c <- cem(treatment = "partisan_elections2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita",  "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "partisan_elections2",  "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat3c$breaks
data2_model3<-data2
data2_model3$match<-mat3c$matched
data2_model3$weight<-mat3c$w
data2_model3<-subset(data2_model3, match=="TRUE")

summary(att(mat3c,taxes_capita~mrp_ideology*partisan_elections2+median_income+percent_black+city_pop+house_value+factor(abb), data=data2))
partisanelections_taxes_pc<-lm(taxes_capita~mrp_ideology*partisan_elections2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model3)

##salestax_share
data2<-subset(data, !is.na(salestax_share))
mat3d.st <- cem(treatment = "partisan_elections2", data = data2[data2$local_salestax_allowed==1,], drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita",  "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "average_payroll_fulltime","partisan_elections2",   "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))

att(mat3d.st,salestax_share~mrp_ideology*partisan_elections2 +median_income+percent_black+city_pop+house_value+factor(abb), data=data2[data2$local_salestax_allowed==1,])
partisanelections_salestaxes<-lm(salestax_share~mrp_ideology*partisan_elections2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model3[data2_model3$local_salestax_allowed==1,])


#####
## initiative (#3)
#####

## Policy
data2<-subset(data, !is.na(policy))
mat4a <- cem(treatment = "initiative2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "welfare_capita", "police_capita","local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "average_payroll_fulltime",  "initiative2", "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat4a$breaks
data2_model4<-data2
data2_model4$match<-mat4a$matched
data2_model4$weight<-mat4a$w
data2_model4<-subset(data2_model4, match=="TRUE")

att(mat4a,policy~mrp_ideology*initiative2 +median_income+city_pop, data=data2)
initiative_policy<-lm(policy~mrp_ideology*initiative2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model4)

## expenditures_capita
data2<-subset(data, !is.na(expenditures_capita))
mat4b <- cem(treatment = "initiative2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "welfare_capita", "police_capita","local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "average_payroll_fulltime",   "initiative2","mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat4b$breaks
data2_model4<-data2
data2_model4$match<-mat4b$matched
data2_model4$weight<-mat4b$w
data2_model4<-subset(data2_model4, match=="TRUE")

att(mat4b,expenditures_capita~mrp_ideology*initiative2 +median_income+city_pop, data=data2)
initiative_exp_pc<-lm(expenditures_capita~mrp_ideology*initiative2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model4)

## taxes_capita
data2<-subset(data, !is.na(taxes_capita))
mat4c <- cem(treatment = "initiative2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "welfare_capita", "police_capita","local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "average_payroll_fulltime",  "initiative2", "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat4c$breaks
data2_model4<-data2
data2_model4$match<-mat4c$matched
data2_model4$weight<-mat4c$w
data2_model4<-subset(data2_model4, match=="TRUE")

att(mat4c,taxes_capita~mrp_ideology*initiative2 +median_income+city_pop, data=data2)
initiative_taxes_pc<-lm(taxes_capita~mrp_ideology*initiative2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model4)

## salestax_share
data2<-subset(data, !is.na(salestax_share))
mat4d.st <- cem(treatment = "initiative2", data = data2[data2$local_salestax_allowed==1,], drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "welfare_capita", "police_capita","local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "average_payroll_fulltime","initiative2",   "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))

att(mat4d.st,salestax_share~mrp_ideology*initiative2 +median_income+city_pop+house_value, data=data2[data2$local_salestax_allowed==1,])
initiative_salestaxes<-lm(salestax_share~mrp_ideology*initiative2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model4[data2_model4$local_salestax_allowed==1,])

#####
## term limits (#4)
#####

## Policy
data2<-subset(data, !is.na(policy))
mat5a <- cem(treatment = "term_limits2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "welfare_capita", "police_capita","local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "average_payroll_fulltime",  "term_limits2",   "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat5a$breaks
data2_model5<-data2
data2_model5$match<-mat5a$matched
data2_model5$weight<-mat5a$w
data2_model5<-subset(data2_model5, match=="TRUE")

att(mat5a,policy~mrp_ideology*term_limits2 +median_income+city_pop, data=data2)
termlimits_policy<-lm(policy~mrp_ideology*term_limits2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model5)

## expenditures_capita
data2<-subset(data, !is.na(expenditures_capita))
mat5b <- cem(treatment = "term_limits2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "welfare_capita", "police_capita","local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "average_payroll_fulltime",   "term_limits2",   "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat5b$breaks
data2_model5<-data2
data2_model5$match<-mat5b$matched
data2_model5$weight<-mat5b$w
data2_model5<-subset(data2_model5, match=="TRUE")

att(mat5b,expenditures_capita~mrp_ideology*term_limits2 +median_income+city_pop, data=data2)
termlimits_exp_pc<-lm(expenditures_capita~mrp_ideology*term_limits2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model5)

## taxes_capita
data2<-subset(data, !is.na(taxes_capita))
mat5c <- cem(treatment = "term_limits2", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "welfare_capita", "police_capita","local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "average_payroll_fulltime",   "term_limits2",   "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat5c$breaks
data2_model5<-data2
data2_model5$match<-mat5c$matched
data2_model5$weight<-mat5c$w
data2_model5<-subset(data2_model5, match=="TRUE")

att(mat5c,taxes_capita~mrp_ideology*term_limits2 +median_income+city_pop, data=data2)
termlimits_taxes_pc<-lm(taxes_capita~mrp_ideology*term_limits2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model5)

## salestax_share
data2<-subset(data, !is.na(salestax_share))
mat5d.st <- cem(treatment = "term_limits2", data = data2[data2$local_salestax_allowed==1,], drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "welfare_capita", "police_capita","local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",  "average_payroll_fulltime", "term_limits2",  "termlength2",   "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))

att(mat5d.st,salestax_share~mrp_ideology*term_limits2 +median_income+city_pop+house_value, data=data2[data2$local_salestax_allowed==1,])
termlimits_salestax<-lm(salestax_share~mrp_ideology*term_limits2 +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model5[data2_model5$local_salestax_allowed==1,])

#####
## Percent at large (#5)
#####

## policy
data2<-subset(data, !is.na(policy))
mat6a <- cem(treatment = "pal_binary", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",     "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black", "perc_atlarge2", "pal_binary"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat6a$breaks

data2_model6<-data2
data2_model6$match<-mat6a$matched
data2_model6$weight<-mat6a$w
data2_model6<-subset(data2_model6, match=="TRUE")

summary(lm(policy~mrp_ideology*pal_binary +percent_black+median_income+city_pop+house_value, data = data2_model6),weights=weight)
pal_policy<-att(mat6a,policy~mrp_ideology*pal_binary +percent_black+median_income+city_pop+house_value+factor(abb), data = data2)
pal_policy<-lm(policy~mrp_ideology*pal_binary +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model6)

## expenditures_capita
data2<-subset(data, !is.na(expenditures_capita))
mat6b <- cem(treatment = "pal_binary", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",     "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black", "perc_atlarge2", "pal_binary"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat6b$breaks

data2_model6<-data2
data2_model6$match<-mat6b$matched
data2_model6$weight<-mat6b$w
data2_model6<-subset(data2_model6, match=="TRUE")

att(mat6b,expenditures_capita~mrp_ideology*pal_binary +percent_black+median_income+city_pop+house_value+factor(abb), data = data2)
pal_exp_pc<-lm(expenditures_capita~mrp_ideology*pal_binary +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model6)

## taxes_capita
data2<-subset(data, !is.na(taxes_capita))
mat6c <- cem(treatment = "pal_binary", data = data2, drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",     "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black", "perc_atlarge2", "pal_binary"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))
mat6c$breaks

data2_model6<-data2
data2_model6$match<-mat6c$matched
data2_model6$weight<-mat6c$w
data2_model6<-subset(data2_model6, match=="TRUE")

att(mat6c,taxes_capita~mrp_ideology*pal_binary +percent_black+median_income+city_pop+house_value+factor(abb), data = data2)
pal_taxes_pc<-lm(taxes_capita~mrp_ideology*pal_binary +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model6)

## salestax_share
data2<-subset(data, !is.na(salestax_share))
mat6d.st <- cem(treatment = "pal_binary", data = data2[data2$local_salestax_allowed==1,], drop = c("city_id","gov_id", "abb", "city.x", "place_fips", "policy", "expenditures_capita", "taxes_capita", "local_salestax_allowed", "salestax_share", "pres_2008", "expenditures_capita",     "mrp_ideology_sd", "mrp_ideology_lower", "mrp_ideology_upper", "percent_black", "perc_atlarge2", "pal_binary"), baseline.group=0,cutpoints = list(house_value=house_cut2,city_pop= pop_cut2,median_income= income_cut2,  mrp_ideology=mrp_cut2))

att(mat6d.st,salestax_share~mrp_ideology*pal_binary +percent_black+median_income+city_pop+house_value+factor(abb), data = data2[data2$local_salestax_allowed==1,])
pal_salestax<-lm(salestax_share~mrp_ideology*pal_binary +median_income+percent_black+city_pop+house_value+factor(abb), weight=weight, data=data2_model6[data2_model6$local_salestax_allowed==1,])

## output the CEM regression models

texreg(list(fog_policy, fog_exp_pc,fog_taxes_pc, fog_salestax))
texreg(list(partisanelections_policy,partisanelections_exp_pc, partisanelections_taxes_pc, partisanelections_salestaxes))
texreg(list(initiative_policy, initiative_exp_pc, initiative_taxes_pc,initiative_salestaxes))
texreg(list(termlimits_policy,termlimits_exp_pc, termlimits_taxes_pc, termlimits_salestax))
texreg(list(pal_policy,pal_exp_pc, pal_taxes_pc, pal_salestax))



###################
###Make CEM Plots
###################

dependentvar <- c("policy","taxes_capita","salestax_share","expenditures_capita")
depvarnames <- c("Policy Scale","Taxes per Capita","Share of Taxes From Sales Taxes","Expenditures per Capita")

matchingplot <- function(dependentvar,depvarnames,datname,match,institution,span=.75,...){
    for (i in dependentvar){
        dat <- get(datname)
        matnames <- c("a","c","d.st","b")
        print(paste(match,matnames[dependentvar==i],sep=""))
        mat <- get(paste(match,matnames[dependentvar==i],sep=""))
        dvar <- dat[,i]
        dat<-subset(get(datname), ! is.na(dvar))
        if(i=="salestax_share"){dat <- dat[dat$local_salestax_allowed==1,]}
        dat$match<-mat$matched
        dat$weight<- mat$w
        dat <- subset(dat, match=="TRUE")
        dvar <- dat[,i]
        
        inst <- dat[[institution]]
        colors <- ifelse(inst==1,"black","grey")
        colors[is.na(colors)] <- "white"
      par(mar = c(5,6,4,2) + 0.1) 
   plot(dvar~dat[["mrp_ideology"]],axes=F,
        ylab="",xlab="City Policy Conservatism",
        ,col=colors,pch=20)
  title(ylab =depvarnames[dependentvar==i], cex.lab = 1,
      line = 3.25)
      
       axis(1,at=c(-1,-.5,0,.5))
       axis(2,las=2)


        temp <- data.frame(x=dat[["mrp_ideology"]][inst==1], y=dvar[inst==1])
        print(dim(temp))
        y.loess <- loess(y ~ x,temp,weight=dat$weight[inst==1],span=span)
        y.predict <- predict(y.loess, data.frame(x=seq(min(temp$x,na.rm=T),max(temp$x,na.rm=T),sd(temp$x,na.rm=T)/5)))
        lines(y.predict~seq(min(temp$x,na.rm=T),max(temp$x,na.rm=T),sd(temp$x,na.rm=T)/5),col="black")
        
        temp <- data.frame(x=dat[["mrp_ideology"]][inst==0], y=dvar[inst==0])
        print(dim(temp))
        y.loess <- loess(y ~ x,temp,weight=dat$weight[inst==0],span=span)
        y.predict <- predict(y.loess, data.frame(x=seq(min(temp$x,na.rm=T),max(temp$x,na.rm=T),sd(temp$x,na.rm=T)/5)))
        lines(y.predict~seq(min(temp$x,na.rm=T),max(temp$x,na.rm=T),sd(temp$x,na.rm=T)/5),col="grey")
    }
}

#############
## Figure 15: Type of Government
#############
pdf("figure15_gov2weightGrey2.pdf")
par(mfrow=c(2,2))
match <- "mat2"
institution <- "fog2"
matchingplot(dependentvar,depvarnames,datname="data",match,institution,span=.85,cex=.75)
dev.off()

#############
## Figure 16: Partisan Elections
#############
pdf("figure16_partisan_elections2weightGrey2.pdf")
par(mfrow=c(2,2))
match <- "mat3"
institution <- "partisan_elections2"
matchingplot(dependentvar,depvarnames,datname="data",match,institution,span=.85)
dev.off()

#############
## Figure 17: Direct Democracy
#############
pdf("figure17_initiative2weightGrey2.pdf")
par(mfrow=c(2,2))
match <- "mat4"
institution <- "initiative2"
matchingplot(dependentvar,depvarnames,datname="data",match,institution,span=1,cex=.75)
dev.off()

#############
## Figure 18: Term Limits
#############
pdf("figure18_term_limits2weightGrey2.pdf")
par(mfrow=c(2,2))
match <- "mat5"
institution <- "term_limits2"
matchingplot(dependentvar,depvarnames,datname="data",match,institution,span=1,cex=.5)
dev.off()

#############
## Figure 19: At-large elections
#############
pdf("figure19_pal_binaryweightGrey2.pdf")
par(mfrow=c(2,2))
match <- "mat6"
institution <- "pal_binary"
matchingplot(dependentvar,depvarnames,datname="data",match,institution,span=1,cex=.5)
dev.off()



